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ABSTRACT 

In  most  sensor  network  applications,  the  vector  containing  the  ob¬ 
servations  gathered  by  the  sensors  lies  in  a  space  of  dimension  equal 
to  the  number  of  nodes,  typically  because  of  observation  noise,  even 
though  the  useful  signal  belongs  to  a  subspace  of  much  smaller  di¬ 
mension.  This  motivates  smoothing  or  rank  reduction.  We  formulate 
a  convex  optimization  problem,  where  we  incorporate  a  fidelity  con¬ 
straint  that  prevents  the  final  smoothed  estimate  from  diverging  too 
far  from  the  observations.  This  leads  to  a  distributed  algorithm  in 
which  nodes  exchange  updates  only  with  neighboring  nodes.  We 
show  that  the  widely  studied  consensus  algorithm  is  indeed  only  a 
very  specific  case  of  our  more  general  formulation.  Finally,  we  study 
the  convergence  rate  and  propose  some  approaches  to  maximize  it. 

Index  Terms —  Distributed  smoothing,  consensus  algorithm, 
sensor  networks. 

1.  INTRODUCTION 

In  this  paper  we  consider  a  wireless  sensor  network  that  monitors  a 
field  of  values,  e.g.,  a  field  of  temperatures  or  the  concentration  of 
a  given  contaminant.  In  most  cases,  the  useful  signal  is  a  smoothed 
function,  as  a  result  of  a  diffusion  process.  However,  typically  the 
set  of  measurements  is  not  at  all  smoothed  because  of  the  observa¬ 
tion  noise.  This  motivates  the  design  of  sensor  networks  that  are  able 
to  apply  some  kind  of  filtering  or  smoothing  to  the  data,  in  order  to 
limit  the  influence  of  noise  on  the  final  decision.  If  possible,  it  would 
be  important  to  implement  these  filtering  operations  in  a  distributed 
way,  i.e.  without  the  presence  of  a  sink  node  that  has  to  collect  all 
the  data  to  apply  the  required  filtering  operation.  This  problem  has 
been  studied  extensively  in  the  case  where  the  useful  signal  is  ho¬ 
mogeneous,  that  is  spatially  constant.  In  such  a  case,  the  so  called 
consensus  algorithms  are  able  to  provide  the  globally  optimal  es¬ 
timate  with  a  network  of  only  locally  interacting  sensors,  see  e.g., 
[1],  [2],  Another  problem  that  has  received  significant  attention  has 
been  the  so  called  distributed  Kalman  filtering  algorithm,  see,  e.g. 
[3],  [4],  where  the  goal  is  to  track  a  non-stationary  process,  typically 
modeled  as  a  Gauss-Markov  field. 

Conversely,  the  goal  of  this  work  is  to  provide  a  distributed 
mechanism  to  perform  a  spatial  smoothing  of  a  stationary,  inhomo¬ 
geneous  field.  Furthermore,  we  incorporate  a  fidelity  constraint  in 
the  optimization  process  to  prevent  the  the  final  estimate  from  di¬ 
verging  too  far  from  the  observations.  In  particular,  we  show  that  if 
the  useful  signal  is  a  continuous  function,  i.e.,  it  can  be  locally  ap¬ 
proximated  by  a  low  order  polynomial  expansion,  then  the  smooth- 
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ing  operation  can  be  implemented  with  a  totally  distributed  network 
where  each  node  is  linked  only  to  a  set  of  neighbors,  with  the  degree 
of  each  node  depending  only  on  the  degree  of  the  approximating 
polynomial.  Then  we  generalize  the  algorithm  to  handle  the  case  of 
discontinuous  functions.  As  expected,  we  will  show  that  the  well 
known  consensus  algorithm  can  be  seen  as  a  particular  case  of  our 
more  general  formulation. 

The  paper  is  organized  as  follows.  In  Section  2,  we  formulate 
the  problem  and  motivate  its  usefulness.  In  Section  3,  we  propose  a 
spatial  smoothing  mechanism  with  no  fidelity  constraint.  In  Section 
4,  we  add  the  fidelity  constraint.  Finally,  in  Section  5,  we  show  how 
to  maximize  the  convergence  rate. 

2.  MOTIVATIONS  AND  PROBLEM  STATEMENT 

We  consider  a  sensor  network  composed  of  N  nodes,  each  equipped 
with  three  basic  components:i)  a  transducer  that  senses  the  physical 
parameter  of  interest;  ii)  a  discrete  dynamical  system  whose  state 
is  initialized  with  the  local  measurements;  iii)  a  radio  interface  that 
transmits  the  state  of  the  dynamical  system  and  receives  the  state 
transmitted  by  the  other  nodes,  thus  ensuring  interaction  among  the 
nodes.  The  observation  of  the  sensor  located  at  (xt,  yi )  is 

9(xi,yi)  =  f(xi,yi)  +w(xi,yi),  i  =  l,...,N  (1) 

where  f(xt,  yi)  is  the  useful  field  and  w(xi,yi)  is  noise.  In  vector 
notation,  we  may  write  g  =  f  +  w,  with  the  three  vectors  hav¬ 
ing  dimension  N.  For  simplicity,  we  assume  that  the  nodes  are 
uniformly  spaced  over  a  2-dimensional  grid  and  that  the  observed 
field  does  not  vary  with  time.  However,  the  extension  to  the  more 
general  case  can  be  done  in  a  straightforward  manner.  Our  goal  is 
to  provide  a  smoothed  version  f(x,  y)  of  the  useful  field  f(x,  y), 
in  order  to  attenuate  as  much  as  possible  the  effect  of  noise,  using 
a  totally  distributed  algorithm.  If  the  useful  signal  is  constant,  its 
value  can  be  obtained,  for  example,  with  the  well  known  average 
consensus  (AC)  algorithm,  which  is  totally  distributed.  The  AC  al¬ 
gorithm  proceeds  as  follows.  Every  node  initializes  its  own  state 
value  Zi (0)  =  g(xt,  yi)  and  then  it  updates  this  value  through  linear 
interactions  with  other  nodes  according  to  the  following  equation 

Zi(k  +  1)  =  Zi(k)  +  s  ^2  an{zi(k)  -  Zi{k)),  (2) 

jeATi 

where  AT%  is  the  set  of  neighbors  of  node  i,  the  coefficients  a,,  are 
real  and  non-negative  (an  =  0),  and  k  is  the  time  or  iteration  in¬ 
dex.  Although  synchronous  updates  are  assumed  here,  these  results 
extend  to  the  asynchronous  case,  and  to  stationary  node  and  link 
failure  models  easily.  Equation  (2)  can  be  recast  in  matrix  form  as 
[1] 

z(k  +  !)  =  (/  —  eL)z(k),  (3) 
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where  L  =  D  —  A  is  the  graph  Laplacian,  A  is  the  adjacency 
matrix,  with  entries  the  coefficients  and  D  is  the  degree  matrix, 
i.e.  a  diagonal  matrix  whose  i-th  entry  is  di  =  YljeM-  a*j-  ^or 
purpose  of  this  work,  the  matrix  A  (and  then  L)  may  be  assumed 
to  be  a  symmetric  matrix.  Let  Aj w(L)  denote  the  largest  eigenvalue 
of  L.  It  is  straightforward  to  show  (see,  e.g.  [1])  that,  if  0  <  e  < 
2/Aj v(L)  and  the  graph  describing  the  network  is  connected  and 
balanced,  then  the  vector  z(k)  converges  to  the  average  consensus, 
i.e., 

lim  z(fc)  =  ^HT9,  (4) 

k — »oo  i\ 

where  1  is  the  vector  with  all  ones. 

This  algorithm  is  totally  decentralized,  but  it  represents  an 
extreme  form  of  smoothing,  as  it  destroys  any  potential  spatial 
variation  in  the  field  of  interest,  which  in  most  cases  is  spatially 
inhomogeneous.  The  goal  of  this  paper  is  to  propose  an  algorithm 
to  smooth  the  observation,  using  a  totally  distributed  approach,  but 
without  forcing  the  final  result  to  converge  to  a  constant  value. 
We  formulate  our  problem  as  an  optimization  problem  with  a  fi¬ 
delity  constraint.  Formally,  we  want  to  find  the  estimated  field 
f  :=  (f(xi,yi), . .  ■ ,  f(xN,  Vn))  as  the  minimizer  of  the  functional 


To  take  into  account  border  effects,  A f  is  the  set  of  indices  for  which 
the  above  differences  can  be  properly  computed;  also  recall  that  we 
have  assumed  a  uniform  2D  grid  for  simplicity  of  exposition.  Since 
(6)  is  a  quadratic  form  on  g,  its  minimum  can  be  reached  using  the 
steepest  descent  method 

z(k  +  1)  =  z(k)  —  eLz(k)  =  Wz(k),  (8) 

with  initialization  z  (0)  =  g,  having  set  W  =  I  —  eL.  It  is  useful  to 
remark  that,  by  construction  L  is  a  positive  semidefinite,  symmetric, 
band  matrix.  Hence,  there  exists  a  unitary  matrix  U  and  a  diagonal 
matrix  A  with  real  elements  such  that: 


where  U2  has  dimension  N  x  L,  with  L  denoting  the  dimension  of 
the  kernel  of  L.  The  columns  of  U2  are  the  vectors  spanning  the 
kernel  of  L.  Let  Ai(L)  and  A;(W)  denote  the  eigenvalues  of  L  and 
W;  we  assume  that  these  eigenvalues  are  ordered  in  non-decreasing 
order.  We  can  always  choose  s,  so  that  the  eigenvalues  of  W  satisfy 
0  <  |Ai(W)|  <  1,  VI.  This  property  is  clearly  achieved  by  setting 


J(f)=^V(f)  +  (l-M)||f-9||2,  (5) 


0  <  £  < 


2 

Ajv(L) 


(10) 


where  V (f)  is  a  convex  function  measuring  the  smoothness  of  f,  the 
term  ||f  —  g||2  measures  the  distance  between  the  final  estimate  and 
the  initial  observation  g,  and  g  is  a  real  parameter  lying  between  0 
and  1.  By  varying  g,  we  vary  the  relative  importance  of  smoothing 
versus  that  of  fidelity  to  the  initial  measurement.  Since  J  (f)  is  a  con¬ 
vex  function  off,  the  minimization  of  J  (f)  admits  a  unique  solution. 
What  is  critical  in  our  approach  is  that  the  minimization  of  (5)  has 
to  be  achieved  using  a  network  with  no  fusion  center,  whose  nodes 
exchange  information  only  with  its  neighbors.  This  is  useful  to  sim¬ 
plify  the  operations  at  each  node  and,  most  important,  to  limit  the 
waste  of  energy  related  to  the  exchange  of  data  among  sensors.  We 
will  consider  first  the  case  of  pure  smoothing  (no  fidelity  constraint). 
Then,  we  will  add  the  fidelity  constraint. 


With  this  choice,  it  is  straightforward  to  verify  that 

AT 

lim  z (k)  =  lim  V'  \\ (W)u,u*z(0)  =  U2U^z(0),  (11) 

k — *oo  k — *00  *  ^ 

i=  1 

where  the  columns  of  U 2  are  exactly  the  vectors  spanning  the  kernel 
of  L.  Hence,  expression  (11)  states  that  the  final  value  coincides 
with  the  projection  of  the  observation  vector  onto  the  nullspace 
of  L.  This  is  an  important  property  that  significantly  enlarges  the 
possibility  of  smoothing  or,  more  generally,  subspace  reduction  of 
the  observation,  by  a  proper  construction  of  the  matrix  L.  Some 
examples  may  be  helpful  to  grasp  the  structural  properties  of  L  and, 
as  a  consequence,  the  final  value. 


3.  SPATIAL  SMOOTHING  WITH  NO  FIDELITY 
CONSTRAINT 


Let  us  consider  the  one-dimensional  case,  for  simplicity  of  notation. 
In  such  a  case,  we  have,  for  example: 


If  the  useful  field  /( x,  y)  is  a  continuous  function  of  the  spatial  co¬ 
ordinates,  according  to  the  Weierstrass’  theorem  it  can  be  approxi¬ 
mated  by  a  two-dimensional  polynomial  of  finite  order  in  the  vari¬ 
ables  x  and  y,  with  an  arbitrarily  small  error.  Let  us  denote  with  K 
the  order  of  the  polynomial  in  both  variables  x  and  y.  In  such  a  case, 
we  want  to  minimize  the  following  function 

V{g)  =  \  £  £  [V<K-mMm)v(-^)]2,  (6) 

i€A/"  m= 0 


a)  K  =  1, 

/  1  -1  0  0 

-1  2-1 


L  = 


0  ^ 

0 

0 


V  0 


-1  2-1 

0-1  1  ) 


where  virn'1  and  Vym)  denote  the  m—th  order  difference  operator 
with  respect  to  the  variables  x  and  y,  respectively.  More  specifically, 
the  operator  is  defined  through  the  following  properties: 

=g(xi,yi);  vi0)g(xi,yi )  =  g(xi,yi); 

Vi  j(*t,  yf)  =  g(xi,yt)  -  g(xi-i,yi); 

Vy 1]g{xi,yt)  =  g(xi,yf)  -  g(xt,  yt-i) 

^(x]g{xi,yi)  =  \y('x~1) g{xi,yf^  ;... 


In  this  case,  the  role  of  L  is  the  same  as  the  Laplacian  matrix  in 
(3),  where  a,j  =  1,  if  \i  —  j\  =  1,  and  0  otherwise.  This  happens 
when  each  node  has  only  two  neighbors  (except  the  border  nodes 
having  only  one  neighbor).  The  matrix  L  has,  in  this  case,  a  null 
eigenvalue  of  multiplicity  one.  Since  each  row  of  L  has  zero  row 
sum,  the  eigenvector  associated  with  the  null  eigenvalue  of  L  is  the 
vector  1  composed  of  all  ones.  Hence,  the  final  result  is  as  in  (4),  as 
in  the  conventional  average  consensus  algorithm. 
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b)  K  =  2, 


/ 

1 

-2 

1 

0  0  0 

0 

-2 

5 

-4 

0 

1 

-4 

6 

0 

0 

6 

-4 

1 

■■  -4 

5 

-2 

V 

0 

0  1 

-2 

1  / 

In  this  case,  the  nullspace  of  L  has  dimensionality  two  and  it  is 
spanned  hy  a  linear  combinations  of  polynomials  of  degree  zero  and 
one.  An  orthonormal  set  is  given,  in  this  case,  by  the  Legendre  poly¬ 
nomials  of  degree  zero  and  one.  Hence  the  final  vector  is  a  straight 
line.  In  general,  it  is  easy  to  check  that,  given  the  structure  of  the 
cost  function  V (g)  in  (6),  the  nullspace  of  L  has  dimensionality  K 
and  it  is  spanned  by  the  polynomials  of  orders  from  0  to  K  —  1. 
Since  any  continuous  function  can  be  approximated  with  an  arbitrar¬ 
ily  small  error,  hy  a  polynomial,  the  above  method  provides  then  a 
distributed  tool  to  approximate  any  continuous  field  of  values.  The 
network  degree  corresponding  to  a  given  value  of  K  is  simply  2 K. 

An  example  of  the  application  of  the  previous  method  to  a  2D 
field  is  reported  in  Fig.  1 .  In  this  case,  the  useful  field  is  a  paraboloid 
and  the  measurement  is  corrupted  hy  zero  mean  white  Gaussian 
noise.  The  observations  are  represented,  in  Fig.  1,  hy  the  circles. 
The  final  state  vector  is  represented  by  the  paraboloid  shown  in  Fig. 
1,  which  is  almost  perfectly  superimposed  on  the  useful  field. 


Fig.  1.  Reconstruction  of  a  2D  field. 


4.  SPATIAL  SMOOTHING  WITH  FIDELITY  CONSTRAINT 

In  principle,  we  can  improve  the  approximation  in  the  above  method 
hy  increasing  the  value  of  K.  However,  the  degree  of  the  network, 
defined  as  the  maximum  number  of  neighbors  of  each  node,  in¬ 
creases  with  K.  This  induces  a  greater  waste  of  energy.  Hence,  it 
is  necessary  to  find  the  right  trade-off  between  energy  consumption 
and  approximation  error.  Furthermore,  if  the  useful  field  presents 
discontinuities,  we  could  be  forced  to  use  very  high  values  of  K  and 


still  have  a  nonnegligible  error.  It  is  then  useful  to  devise  some  vari¬ 
ants  of  the  previous  method  that  still  allow  us  to  have  a  limited  value 
of  the  network  degree.  To  this  end,  we  use  now  as  cost  function  the 
expression  (5),  with  V ( g )  given  in  (6),  with  g  £  (0, 1).  In  this  case, 
the  steepest  descent  method  leads  to 

z(k  +  1)  =  [(1  —  e(l  —  g))I  —  £gh]z(k)  +  e(l  —  A) z(0),  (12) 

with  z( 0)  =  g.  Introducing  the  matrix  W  =  (1  —  e(l  —  g))I  — 
£gL,  we  can  guarantee  that  the  eigenvalues  Ai(W)  of  W  are  strictly 
between  —  1  and  1,  by  setting 

(21  2 

0  <  e  <  min  <  - - - >  =  - — - - - .  (13) 

i  \  \i(L)  -J-  1  —  g  J  A n(L)  +  1  —  g 

With  a  few  simple  algebraic  manipulations,  we  can  rewrite  z(k )  as 

z(k)  =  Wfcz(0)+e(l-A*)El^Wfc*(0) 

=  Wkz(0)  +  e(l  -  p)(I  -  W)-1^  -  W k)z(k). 

(14) 

Choosing  e  according  to  (13),  the  state  vector  converges  to 

lim  z(k)  =  (l+  -1 —  l)  g.  (15) 

fc— >oo  y  I  —  fj,  J 

Depending  on  the  value  of  g,  we  may  give  different  relative  impor¬ 
tance  to  smoothing  or  fidelity  to  the  original  observation.  In  the  ex¬ 
treme  case  of  g  =  0,  the  network  does  not  apply  any  smoothing,  i.e. 
limjc^oo  z(k)  =  g,  whereas,  at  the  other  extreme,  when  g  =  1,  the 
final  value  coincides  with  the  projection  of  the  observation  onto  the 
nullspace  of  the  Laplacian  matrix,  as  proved  in  the  previous  section. 

A  numerical  example  is  reported  in  Fig.  2,  relative  to  a  one¬ 
dimensional  network  located  over  a  straight  line.  The  observed  sig¬ 
nal  in  this  case  is  a  sinusoid  (dashed  line)  and  the  observation  (dots) 
is  corrupted  by  white  Gaussian  noise.  The  SNR  is  5  dB.  Smoothing 
has  been  performed  using  the  simple  algorithm  (8),  with  I\  =  3. 
In  this  case,  with  g  =  .1,  the  method  projects  the  observed  vec¬ 
tor  onto  the  space  spanned  by  second  order  polynomials.  Since  the 
observation  is  a  noisy  sinusoid,  the  final  result  (dash-dotted  line)  is 
not  very  good.  However,  as  soon  as  g  is  slightly  less  than  one,  the 
method  is  forced  to  take  into  account  the  fidelity  to  the  observation, 
and  the  final  result  is  much  better  than  in  the  previous  case.  Using 
g  =  0.9999,  for  example,  the  result  of  the  smoothing  operation  is 
represented  by  the  solid  line  and  we  can  see  that  the  approximation 
is  now  pretty  good. 

Even  if  this  is  only  a  simple  example,  Fig.  2  suggests  that  the 
choice  of  g  can  have  a  strong  impact  on  the  smoothing  operation. 
To  quantify  the  final  distortion,  we  can  compute  the  mean  square 
error,  averaged  over  the  noise  realizations.  Introducing  the  matrix 

P(g)  —  ^1  +  L'j  ,  the  final  MSE  is 

MSE(g)  =  ||(P(M)  -  I)f\\2  +  a2ntr(P(g)P(g)T).  (16) 

In  the  case  of  a  sinusoidal  function,  this  function,  normalized  to 
||/||2,  is  reported  in  Fig.  3,  for  different  values  of  g  and  cr2.  As 
expected,  there  is  an  optimal  value  of  g  that  depends  on  the  noise 
level:  When  there  is  no  noise,  it  is  better  to  apply  no  smoothing  at 
all,  and  thus  the  best  value  of  g  is  zero;  conversely,  as  the  noise 
increases,  it  is  better  to  use  values  of  g  closer  and  closer  to  one. 
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Fig.  2.  Reconstruction  of  a  noisy  sinusoid. 


Fig.  3.  MSE  as  a  function  of  n  and  an. 


5.  CONVERGENCE  RATE 

A  critical  issue  in  most  distributed  processing  algorithms  to  be  im¬ 
plemented  in  decentralized  sensor  networks  is  the  convergence  rate. 
It  is  clear  in  fact  that  a  slow  convergence  rate  implies  a  waste  of  time 
and  energy.  And  energy  is  an  especially  critical  and  scarce  resource. 

The  distributed  projection  method  illustrated  in  Section  3  re¬ 
quires  only  that  the  nullspace  of  the  matrix  L  to  be  used  for  the 
state  update  in  (8)  be  composed  of  the  right  eigenvectors  (polyno¬ 
mials),  hut  it  does  not  impose  any  condition  on  the  eigenvalues  of 
( L ).  Hence,  we  can  optimize  the  choice  of  the  eigenvalues  in  order 
to  maximize  the  convergence  rate. 

From  (1 1),  the  convergence  is  exponential  and  the  convergence 
rate  is  given  by  the  slowest  mode,  i.e.  the  highest  eigenvalue  of 
W  strictly  less  than  1.  This  value  is  Aat_i(W)  =  1  —  5X2 {L). 
Requiring  that  the  corresponding  mode  decays  from  1  to  a,  with 
a  <  1,  we  get  a  convergence  time 

T  —  _ _  (17) 

c  log(l-eA2(L))-  (  ’ 

From  (17),  it  is  clear  that  the  minimization  of  the  convergence  time 
corresponds  to  minimizing  the  second  smallest  eigenvalue  of  L.  To 
make  a  fair  comparison  with  the  method  presented  in  Section  3, 
with  the  matrix  L  resulting  from  the  computation  of  the  gradient 
of  V(g)  in  (6),  we  impose  a  constraint  on  the  trace  of  L  and  look 


for  the  optimal  distribution  of  eigenvalues  that  provides  the  largest 
A2  (L) .  The  result  is  very  simple:  the  optimal  set  of  eigenvalues  cor¬ 
responds  to  having  all  the  eigenvalues  equal  to  each  other  and  equal 
to  tr(i)/rank(£).  We  denote  by  L'  the  matrix  having  the  same 
eigenvectors  as  L,  same  nullspace,  hut  all  the  nonnull  eigenvalues 
equal  to  each  other.  The  price  to  he  paid  for  the  use  of  this  matrix 


Number  of  nodes 


Fig.  4.  MSE  as  a  function  of  fi  and  ern. 

L'  is  a  sort  of  spillover  effect:  Whereas  the  matrix  L  resulting  from 
taking  the  gradient  of  V ( g )  in  (6)  is  a  hand  matrix,  the  matrix  L'  is 
not  necessarily  a  band  matrix.  This  means  that  every  node  in  the  net¬ 
work  has,  in  general,  a  higher  degree  and  this  implies  a  higher  trans¬ 
mit  power  for  each  node.  Setting  L  =  D  —  A,  as  in  the  conventional 
Laplacian,  we  studied  the  spillover  ratio  r(K)  —  Yhx  a%j /  Yb  ai.i > 
where  B  is  the  set  of  indices  falling  in  the  same  band  as  the  initial 
matrix  L ,  whereas  X  indicates  the  set  of  indices  falling  outside  the 
initial  band.  In  Fig.  4,  we  report  the  ratio  of  the  smallest  eigenvalues 
\2{L’) / \2{L)  (solid  lines),  denoting  the  increase  of  convergence 
rate  resulting  from  forcing  all  the  eigenvalues  to  be  equal  to  each 
other,  for  K  =  1,  2,  and  3.  In  the  same  figure,  we  plot  the  spill-over 
ratio  r(K),  also  for  K  =  1,  2,  and  3,  as  a  function  of  the  number 
of  nodes  (dashed  lines).  We  can  notice  how  a  proper  allocation  of 
eigenvalues  produces  a  substantial  increase  of  the  convergence  rate, 
more  and  more  evident  as  the  number  of  nodes  increases.  At  the 
same  time,  this  determines  a  higher  degree  for  each  node,  as  evi¬ 
denced  by  the  dashed  line  curves,  which  implies  a  higher  required 
transmit  power  at  each  node.  An  interesting,  still  open  question,  is 
then  how  to  find  out  the  optimal  eigenvalue  distribution  that  mini¬ 
mizes  the  overall  energy  required  for  achieving  the  final  smoothed 
signal,  taking  into  account  both  convergence  rate  and  local  transmit 
power. 
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